library(RaceID)
library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)
Attaching SeuratObject
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
Loading required package: RColorBrewer
# Seurat files
so_file <- "data/object/seurat.rds"
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so <- readRDS(so_file)
so <- subset(so, subset = nFeature_RNA <= 3000 & nCount_RNA <= 15000)
Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”
so_1 <- subset(so, subset = treatment == "NaCl" & tissue == "Progenitor" & rna_snn_res.0.8 != 18)
so_2 <- subset(so, subset = treatment == "CpG" & tissue == "Progenitor" & rna_snn_res.0.8 != 18)
# library(Seurat)
# Idents(so_1) <- "sample_name"
# deg_1 <- FindMarkers(so_1, ident.1 = "Progenitor_NaCl_Rep1", ident.2 = "Progenitor_NaCl_Rep2")
# Idents(so_2) <- "sample_name"
# deg_2 <- FindMarkers(so_2, ident.1 = "Progenitor_CpG_Rep1", ident.2 = "Progenitor_CpG_Rep2")
sc_1 <- SCseq(expdata = GetAssayData(so_1, slot = "counts", assay = "RNA"))
sc_2 <- SCseq(expdata = GetAssayData(so_2, slot = "counts", assay = "RNA"))
Filterdata can be run with LBatch and bmode to remove genes associated with batch effects. VarID workflow also allows regression of batch effect in pruneKnn with regNB. in that case the regNB parameter of the compNoise function needs to be TRUE as well. Additional variables can be regressed out using the regVar variable with a data.frame. The rownames correspond to cell id and columns to values e.g. cell cycle.
# Filter
sc_1 <- filterdata(sc_1, mintotal = 2500, minexpr = 5, minnumber = 5, CGenes = rownames(so_1)[grep("^(mt)", rownames(so_1))], FGenes = NULL)
sc_2 <- filterdata(sc_2, mintotal = 2500, minexpr = 5, minnumber = 5, CGenes = rownames(so_2)[grep("^(mt)", rownames(so_2))], FGenes = NULL)
expdata_1 <- getExpData(sc_1)
expdata_2 <- getExpData(sc_2)
# Batch
batch_1 <- so_1$sample_name
names(batch_1) <- colnames(so_1)
batch_1 <- batch_1[names(batch_1) %in% colnames(expdata_1)]
batch_2 <- so_2$sample_name
names(batch_2) <- colnames(so_2)
batch_2 <- batch_2[names(batch_2) %in% colnames(expdata_2)]
res_1 <- pruneKnn(expdata_1, large = TRUE, regNB = TRUE, batch = batch_1, knn = 10, alpha = NULL, no_cores = NULL, seed = 12345)
cl_1 <- graphCluster(res_1, pvalue = 0.01)
probs_1 <- transitionProbs(res_1, cl_1)
res_2 <- pruneKnn(expdata_2, large = TRUE, regNB = TRUE, batch = batch_2, knn = 10, alpha = NULL, no_cores = NULL, seed = 12345)
cl_2 <- graphCluster(res_2, pvalue = 0.01)
probs_2 <- transitionProbs(res_2, cl_2)
noise_1 <- compNoise(expdata_1, res_1, regNB = TRUE, batch = batch_1, pvalue = 0.01, no_cores = NULL)
noise_2 <- compNoise(expdata_2, res_2, regNB = TRUE, batch = batch_2, pvalue = 0.01, no_cores = NULL)
sc_1 <- updateSC(sc_1, res = res_1, cl = cl_1, noise = noise_1, flo = .1)
sc_2 <- updateSC(sc_2, res = res_2, cl = cl_2, noise = noise_2, flo = .1)
saveRDS(sc_1, "data/object/components/raceid_progenitor_mt_nacl.rds")
saveRDS(sc_2, "data/object/components/raceid_progenitor_mt_cpg.rds")
saveRDS(res_1, "data/object/components/raceid_progenitor_mt_nacl_res.rds")
saveRDS(res_2, "data/object/components/raceid_progenitor_mt_cpg_res.rds")
saveRDS(cl_1, "data/object/components/raceid_progenitor_mt_nacl_cl.rds")
saveRDS(cl_2, "data/object/components/raceid_progenitor_mt_cpg_cl.rds")
saveRDS(probs_1, "data/object/components/raceid_progenitor_mt_nacl_probs.rds")
saveRDS(probs_2, "data/object/components/raceid_progenitor_mt_cpg_probs.rds")
# sc_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl.rds")
# sc_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg.rds")
# res_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_res.rds")
# res_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_res.rds")
# cl_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_cl.rds")
# cl_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_cl.rds")
# probs_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_probs.rds")
# probs_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_probs.rds")
sc_1 <- compumap(sc_1, spread = 3)
sc_2 <- compumap(sc_2, spread = 3)
options(repr.plot.width = 5, repr.plot.height = 5)
umap_1 <- plotmap(sc_1, um = TRUE)
umap_2 <- plotmap(sc_2, um = TRUE)
## map of transition probabilities
options(repr.plot.width = 5, repr.plot.height = 5)
plotTrProbs(sc_1, probs_1, tp = .5, prthr = 0.01, cthr = 0, um = TRUE)
plotTrProbs(sc_2, probs_2, tp = .5, prthr = 0.01, cthr = 0, um = TRUE)
so_1 <- subset(so_1, cells = rownames(sc_1@umap))
so_2 <- subset(so_2, cells = rownames(sc_2@umap))
all(colnames(so_1) == rownames(sc_1@umap))
all(colnames(so_2) == rownames(sc_2@umap))
so_1$cluster <- sc_1@cluster$kpart
so_2$cluster <- sc_2@cluster$kpart
so_1@reductions$varid_umap <- CreateDimReducObject(embeddings = as.matrix(sc_1@umap), assay = "RNA", key = "varid_umap")
so_2@reductions$varid_umap <- CreateDimReducObject(embeddings = as.matrix(sc_2@umap), assay = "RNA", key = "varid_umap")
Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from varid_umap to varidumap_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to varidumap_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from varid_umap to varidumap_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to varidumap_”
reduction <- "varid_umap"
options(repr.plot.width = 20, repr.plot.height = 10)
dplot_1 <- DimPlot(so_1, reduction = reduction, group.by = "sample_name", label = FALSE) &
ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_2 <- DimPlot(so_1, reduction = reduction, group.by = "cluster", label = TRUE) &
ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "none") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so_1, reduction = reduction, group.by = "main_labels", label = FALSE) &
ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so_1, reduction = reduction, group.by = "fine_labels", label = FALSE) &
ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so_2, reduction = reduction, group.by = "sample_name", label = FALSE) &
ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_6 <- DimPlot(so_2, reduction = reduction, group.by = "cluster", label = TRUE) &
ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "none") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_7 <- DimPlot(so_2, reduction = reduction, group.by = "main_labels", label = FALSE) &
ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_8 <- DimPlot(so_2, reduction = reduction, group.by = "fine_labels", label = FALSE) &
ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + dplot_7 + dplot_8 + plot_layout(ncol = 4)
dplot
fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "fine_delta_score") & theme(aspect.ratio = 1)
fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "fine_delta_score") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot
fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "pRp_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "pRp_RNA") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot
fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "msG0_RNA1") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "msInt_RNA1") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "msCC_diff_RNA") & theme(aspect.ratio = 1)
fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "msG0_RNA1") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "msInt_RNA1") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "msCC_diff_RNA") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot
cluster_main_labels <- ggplot(so_1@meta.data, aes(x = factor(cluster), fill = factor(main_labels, levels = names(color$main_labels))))+
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))
all_main_labels <- ggplot(so_1@meta.data, aes(x = "NaCl", fill = factor(main_labels, levels = names(color$main_labels))))+
geom_bar(stat = "count", position = "fill", width = 1) +
scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
ggtitle("Sort frequency progenitor NaCl") + xlab("") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))
cluster_fine_labels <- ggplot(so_1@meta.data, aes(x = factor(cluster), fill = factor(fine_labels, levels = names(color$fine_labels)))) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))
all_fine_labels <- ggplot(so_1@meta.data, aes(x = "NaCl", fill = factor(fine_labels, levels = names(color$fine_labels)))) +
geom_bar(stat = "count", position = "fill", width = 1) +
scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))
options(repr.plot.width = 20, repr.plot.height = 5)
cluster_main_labels + all_main_labels + cluster_fine_labels + all_fine_labels + plot_layout(width = c(1, 1/11, 1, 1/11), ncol = 4, guides = 'collect') & theme(legend.position = 'bottom')
cluster_main_labels <- ggplot(so_2@meta.data, aes(x = factor(cluster), fill = factor(main_labels, levels = names(color$main_labels))))+
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))
all_main_labels <- ggplot(so_2@meta.data, aes(x = "NaCl", fill = factor(main_labels, levels = names(color$main_labels))))+
geom_bar(stat = "count", position = "fill", width = 1) +
scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
ggtitle("Sort frequency progenitor NaCl") + xlab("") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))
cluster_fine_labels <- ggplot(so_2@meta.data, aes(x = factor(cluster), fill = factor(fine_labels, levels = names(color$fine_labels)))) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))
all_fine_labels <- ggplot(so_2@meta.data, aes(x = "NaCl", fill = factor(fine_labels, levels = names(color$fine_labels)))) +
geom_bar(stat = "count", position = "fill", width = 1) +
scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom") +
guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))
options(repr.plot.width = 20, repr.plot.height = 5)
cluster_main_labels + all_main_labels + cluster_fine_labels + all_fine_labels + plot_layout(width = c(1, 1/11, 1, 1/11), ncol = 4, guides = 'collect') & theme(legend.position = 'bottom')